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Abstract Gravitational-wave detections are enabling measurements of the 
rate of coalescences of binaries composed of two compact objects — neutron 
stars and/or black holes. The coalescence rate of binaries containing neutron 
stars is further constrained by electromagnetic observations, including Galac- 
tic radio binary pulsars and short gamma-ray bursts. Meanwhile, increasingly 
sophisticated models of compact objects merging through a variety of evolu- 
tionary channels produce a range of theoretically predicted rates. Rapid im- 
provements in instrument sensitivity, along with plans for new and improved 
surveys, make this an opportune time to summarise the existing observational 
and theoretical knowledge of compact-binary coalescence rates. 
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1 Introduction 


It has been more than a decade since Abadie et al. (2010) reviewed the lit- 
erature on compact-object merger rate predictions. At that time, a handful 
of ultimately merging double neutron star (NS) systems in the Milky Way 
Galaxy were known from radio pulsar observations, starting with the discov- 
ery of Hulse and Taylor (1975). It was believed that short gamma-ray bursts 
(SGRBs) were associated with mergers of two NSs. But no actual compact 
object mergers were definitively observed. Therefore, most of the rate predic- 
tions, particularly for mergers involving black holes (BHs), were based purely 
on theoretical models (see, e.g., Mandel and O’Shaughnessy 2010 for a review 
of the status of models at that time). 

The intervening decade has completely changed this landscape. On 14 
September, 2015, the advanced Laser Interferometer Gravitational-wave Ob- 
servatory (LIGO, Aasi et al. 2015) detected the first chirp of gravitational 
waves from a binary BH merger, GW150914 (Abbott et al., 2016b). Since then, 
nine more binary BH coalescences were observed during the first and second 
observing runs of advanced LIGO and Virgo detectors (Abbott et al., 2019b). 
The data release from the third observing run brought the total number of con- 
fidently detected binary BH coalescences to approximately 70 (Abbott et al., 
2021d). The first detection of a binary neutron star merger, GW170817 (Ab- 
bott et al., 2017b), was accompanied by a short GRB (Abbott et al., 2017a) 
and a kilonova (Abbott et al., 2017c), firmly establishing the connection be- 
tween these phenomena. A second confident binary neutron star detection, 
GW190425 (Abbott et al., 2020a), followed, along with detections of two NS- 
BH mergers, GW200105 and GW200115 (Abbott et al., 2021a). 

At the same time, significant progress has been made in the theoretical 
modelling of the sources of compact binary coalescences. This has included 
more detailed exploration of the classical channel of isolated binary evolution, 
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typically through the common-envelope phase (van den Heuvel, 1976; Smarr 
and Blandford, 1976; Tutukov and Yungelson, 1993) (or possibly through sta- 
ble mass transfer alone, van den Heuvel et al. 2017; Inayoshi et al. 2017), along 
with significant new developments in the investigation of alternative channels, 
including dynamical interaction in dense stellar environments (Sigurdsson and 
Hernquist, 1993; Kulkarni et al., 1993; Portegies Zwart and McMillan, 2000); 
mergers in binaries of rapidly rotating stars undergoing chemically homoge- 
neous evolution (Mandel and de Mink, 2016; Marchant et al., 2016); mergers 
facilitated by the Lidov-Kozai resonance (Lidov, 1962; Kozai, 1962) in hier- 
archical triple systems; and even mergers of primordial black holes of cosmo- 
logical rather than astrophysical origin (Bird et al., 2016; Ali-Haimoud et al., 
2017). 


The coupled advances in observations and theory are making it possible to 
quantitatively compare rate predictions against observations for the first time 
(e.g., Abbott et al., 2019a, 2021e; Belczynski et al., 2018a; Wysocki et al., 
2018; Neijssel et al., 2019; Belczynski et al., 2020; Farmer et al., 2020; Bav- 
era et al., 2021). Beyond this, they provide critical input into understanding 
the chemical enrichment of the Universe, particularly with elements created 
through r-process nucleosynthesis likely associated with mergers involving NSs 
with their neutron-rich material (Kasen et al., 2017; Cóté et al., 2018; Metzger, 
2019; Kobayashi et al., 2020). The merger rates determine prospects for future 
Earth-based and space-borne gravitational-wave detectors, and will impact de- 
tector design (Amaro-Seoane et al., 2017; Adhikari et al., 2019; Reitze et al., 
2019). Merger rates feed into constraints or predictions on other observables, 
ranging from kilonovae associated with binary NS mergers (e.g., Kasliwal et al., 
2020; Mochkovitch et al., 2021) to gravitational-wave stochastic backgrounds 
(Abbott et al., 2016a) to, possibly, fast radio bursts (e.g., Zhang, 2020). This, 
then, is an opportune time for a review of the current state of the observations 
and predictions of compact object merger rates — and, given the expectation 
of continuing rapid development in this field, it particularly calls for a Living 
Review. 


'This review is organised as follows. We provide an executive summary of the 
observed and theoretically predicted compact-binary merger rates in Sect. 2. 
We then provide more detailed information and discussion of the observations 
in Sect. 3 and of the models in Sect. 4 before concluding in Sect. 5. The collated 
data and code to reproduce all tables and figures in this review are publicly 
available through Broekgaarden and Mandel (2022)!, where we also provide 
supplementary information on how we obtained the data in the tables and 
figures. 


1 The current paper is based on version 6 of the data. 
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2 Executive summary 


We summarise the coalescence rates in three tables below: NS-NS binaries in 
Table 1, NS-BH binaries (where we do not distinguish whether the BH or NS 
formed first) in Table 2, and BH-BH binaries in table 3. 


Coalescence rates are, in general, functions of redshift; we quote current 
local rates at redshift z = 0 per unit source time per unit comoving volume 
in units of Gpc^? yr-!, but caution that these could be up to an order of 
magnitude larger at higher redshifts. Where initially stated in different units, 
we convert these, using, as appropriate, factors of 1.7 x 101° solar blue-light 
luminosities per Milky Way equivalent galaxy (MWEG), 1.17 x 107? MWEG 
per Mpc? (Kopparapu et al., 2008), a globular cluster space density of 2.9 
per Mpc? at z = 0 (Portegies Zwart and McMillan, 2000) and a local su- 
pernova rate of 1.06 x 10? Gpc ? yr-! (Taylor et al., 2014). Of course, such 
simple re-scalings do not account for the dependence of merger rates on the 
star formation history and metallicity (de Freitas Pacheco et al., 2006; Bel- 
czynski et al., 2010; Dvorkin et al., 2016; Neijssel et al., 2019; Chruslinska and 
Nelemans, 2019; Mapelli, 2021). 


In general, we follow original papers in stating the rates; a discussion fol- 
lows in the next two sections. Where uncertainties are stated, we try to indi- 
cate them with + subscripts and superscripts; however, these are not always 
available, and may not be attributable to a specific confidence interval when 
available. Sometimes, we state a range in square brackets without a central 
value. 


The vast literature on compact-object merger rate predictions dating back 
to at least 1979 (Clark et al., 1979) makes a complete historical review implau- 
sible. Therefore, we focus on the latest contributions from each group, except 
where contributions use significantly different methodology or examine the 
impact of different assumptions. We generally eschew papers published before 
2010 as these have typically been superseded by observational and theoretical 
advances. 


We combine rates into several broad categories for ease of viewing, but 
multiple different channels with a broad range of physics may be at play within 
any one category. For example, the category of “nuclear star clusters” may 
include dynamical captures in nuclear clusters without a massive black hole, 
mergers in hierarchical triples involving a massive black hole, or mergers aided 
by an accretion disk in an active galactic nucleus. We do not include some 
very rare channels with predicted rates below the range shown in the figures. 

We do not quote predictions for the number of detections per year for a 
given instrument because such predictions must take into account the variation 
in both the merger rate and the mass distribution of merging binaries with 
redshift for detectors sensitive to cosmological distances, and this information 
is frequently not readily available in the literature. Derivations of how to do 
such calculations are provided by, e.g., Belczynski et al. (2014) and de Mink 
and Mandel (2016); both analytical fits (e.g., Fishbach and Holz 2017) and 
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Table 1: Summary of local NS-NS coalescence rates, part 1: observations (Sup- 
plementary material and Zenodo) 


Source Rate Reference 
[Gpe yr-1] 
Observations: 
Gravitational-wave observations, PDB (ind) 44* 97 GWTC-3; Abbott et al. (2021e) 
Gravitational-wave observations, MS 680771079 GW'TC-3; Abbott et al. (2021e) 
Gravitational-wave observations, BGP 98.0+260;0 GWTC-3; Abbott et al. (2021e) 
Short GRBs [5, 1800] Coward et al. (2012) 
Short GRBs [500, 1500] Petrillo et al. (2013) 
Short GRBs 27071550 Fong et al. (2015) 
Short GRBs, based on GW170817 3527850 Della Valle et al. (2018) 
Short GRBs 1109* 1222 Jin et al. (2018) 
Short GRBs, based on GW170817 1907760 Zhang et al. (2018) 
Short GRBs, based on GW170817, SWIFT  160*700 Dichiara et al. (2020) 
Short GRBs, ~ 150 Zevin et al. (2022) 
Kilonovae lower limit > 8.1 Jin et al. (2016) 
Kilonovae, DES, upper limit < 24000 Doctor et al. (2017) 
Kilonovae, PTF, upper limit < 800 Kasliwal et al. (2017) 
Kilonovae, ATLAS, upper limit < 30000 Smartt et al. (2017) 
Kilonovae, DLT40, upper limit < 99000 Yang et al. (2017) 
Kilonovae, ZTF, upper limit < 900 Andreoni et al. (2021) 
Galactic pulsar binaries x 83072170 O'Shaughnessy and Kim (2010) 
Galactic pulsar binaries 2501330 Kim et al. (2015) 
Galactic pulsar binaries scot, Pol et al. (2020) 
Galactic pulsar binaries 370/7156 Grunthal et al. (2021) 


numerical codes (e.g., Gerosa 2017; Team COMPAS: Riley et al. 2022) are 
available. 


3 Observed rates 


In this section, we discuss the latest observational constraints on compact- 
object coalescence rates. There are direct observations of compact binary merg- 
ers, either through gravitational waves or through electromagnetic transients 
accompanying such mergers, particularly short gamma-ray bursts. Alterna- 
tively, there are observations of binaries that are not yet merging, but are 
expected to do so as compact objects in the future, such as Galactic double 
neutron stars and perhaps some high-mass X-ray binaries. 

In all cases the key challenge in accurately inferring rates from limited 
observations comes from significant selection effects, which may be particularly 
difficult to quantify when the data set is inhomogeneous, and/or when the 
distribution of the intrinsic loudness of the signals is not known. Moreover, 
in some cases there can be uncertainty in the interpretation of the signals. 
Meanwhile, for systems that are not yet merging, there may be uncertainties 
in their future evolution. We discuss these challenges below. 
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Table 1: Summary of local NS-NS coalescence rates, part 2: model predictions 


Source Rate Reference 

Gpc-? yr- 1] 
Models: 
Isolated binary population synthesis, StarTrack 30, 1700] O'Shaughnessy et al. (2010) 
Isolated binary population synthesis, Scenario Machine 1050, 3860] Lipunov and Pruzhinskaya (2014) 
Isolated binary population synthesis, Brussels code < 1.3, 1800] Mennekens and Vanbeveren (2014) 
Isolated binary population synthesis, StarTrack 30, 540 de Mink and Belczynski (2015) 
Isolated binary population synthesis, StarTrack 52, 162 Dominik et al. (2015) 
Isolated binary population synthesis, BSE 240, 1800] Ablimit and Maeda (2018) 
Isolated binary population synthesis, StarTrack 8,50] Belczynski et al. (2018a) 
Isolated binary population synthesis, Star' Track 1.5,631] Chruslinska et al. (2018) 
Isolated binary population synthesis, MOBSE 10,510 Giacobbo and Mapelli (2018) 
Isolated binary population synthesis, Star' Track 24, 68] Klencki et al. (2018) 
Isolated binary population synthesis, COMBINE 2.7,159] Kruckow et al. (2018) 
Isolated binary population synthesis, MOBSE 19,591 Mapelli and Giacobbo (2018) 
Isolated binary population synthesis, COMPAS 61.5, 362] Vigna-Gómez et al. (2018) 
Isolated binary population synthesis, MOBSE 238 Artale et al. (2019) 
Isolated binary population synthesis, MOBSE 12, 400 Baibhav et al. (2019) 
Isolated binary population synthesis, SEVN 70 Boco et al. (2019) 
Isolated binary population synthesis, Star Track 48, 885 Chruslinska et al. (2019) 
Isolated binary population synthesis, BPASS 339, 2178] Eldridge et al. (2019) 
Isolated binary population synthesis, COMPAS 20, 245 Neijssel et al. (2019) 
Isolated binary population synthesis, Star'Track 49.3, 524] Belczynski et al. (2020) 
Isolated binary population synthesis, MOBSE 20, 640 Giacobbo and Mapelli (2020) 
Isolated binary population synthesis, MOBSE 283*97 Santoliquido et al. (2020) 
Isolated binary population synthesis, BPASS 394, 3190] Tang et al. (2020) 
Isolated binary population synthesis, COSMIC 600, 8900] Zevin et al. (2020) 
Isolated binary population synthesis, BSE, 0.4, 1404] Chu et al. (2021) 
Isolated binary population synthesis, BPASS 43, 745] Ghodla et al. (2021) 
Isolated binary population synthesis, Star Track 148, 322 Olejak et al. (2021) 
Isolated binary population synthesis, MOBSE 4.3, 1036.8] Santoliquido et al. (2021) 
Isolated binary population synthesis, BPASS, 27 Briel et al. (2022b) 
Isolated binary population synthesis, COMPAS 0.32, 330] Broekgaarden et al. (2021, 2022) 
Isolated binary population synthesis, Star Track 116, 155 Olejak et al. (2022) 
Hierarchical triples, SecularMultiple 164, 3793] Hamers and Thompson (2019) 
Hierarchical quadruples, MSE 0.8, 30.2 Vynatheya and Hamers (2021) 
Globular cluster dynamics 30 Lee et al. (2010) 
Globular cluster dynamics 0.32, 3.2 Bae et al. (2014) 
Globular cluster dynamics 121 Samsing et al. (2014) 
Globular cluster dynamics, MOCCA 0.02, 0.5 Belczynski et al. (20182) 
Globular cluster dynamics, CMC 0.009, < 25.5] Ye et al. (2020) 
Nuclear star cluster dynamics with SMBH 0.004, 1.4] Antonini and Perets (2012) 
Nuclear star cluster dynamics, with SMBH S 0.02 Petrovich and Antonini (2017) 
Nuclear star cluster dynamics 0.007, 0.1] Belczynski et al. (20182) 
Nuclear star cluster dynamics, with SMBH 0.06, 0.1] Fragione et al. (2019) 
Nuclear star cluster dynamics, with SMBH S 400 McKernan et al. (2020) 
Nuclear star cluster dynamics, with SMBH = [0.15, 0.3] Wang et al. (2021) 
Young star clusters 0.03, 0.15] Ziosi et al. (2014) 
Young/Open star clusters, Nbody7 0.01 — 0.1] Fragione and Banerjee (2020) 


Young star clusters, MOBSE 


+59 
15139 


Santoliquido et al. (2020) 
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Fig. 1: NS-NS merger rates from Table 1. Black lines denote the interval bound- 
aries. Triangles mark upper and lower limits. Colored circles mark inferred or 
simulated values, where multiple values are given in a study. Black circles mark 
the median or mean of a confidence interval, where one is provided. (GitHub) 
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Table 2: Summary of local NS-BH coalescence rates (Supplementary material 


and Zenodo) 


Source Rate Reference 

[Gpc ^? yr^!] 
Observations: 
Gravitational-wave observations, PDB (ind) qm GWTC-3; Abbott et al. (2021e) 
Gravitational-wave observations, MS sgt. GW'TC-3; Abbott et al. (2021e) 
Gravitational-wave observations, BGP 320019 GW'TC-3; Abbott et al. (2021e) 
Galactic pulsar binaries < 1800 Pol et al. (2021) 
Models: 
Isolated binary population synthesis, StarTrack 10, 280] O'Shaughnessy et al. (2010) 
Isolated binary population synthesis, Brussels code 0.06, 800] Mennekens and Vanbeveren (2014) 
Isolated binary population synthesis, Star Track 9,115 de Mink and Belczynski (2015) 
Isolated binary population synthesis, Star' Track 0.04, 20] Dominik et al. (2015) 
Isolated binary population synthesis, BSE 1,160 Ablimit and Maeda (2018) 
Isolated binary population synthesis, MOBSE 5, 180 Giacobbo and Mapelli (2018) 
Isolated binary population synthesis, Star' Track 13.3, 26.7] Klencki et al. (2018) 
Isolated binary population synthesis, COMBINE 2,53] Kruckow et al. (2018) 
Isolated binary population synthesis, MOBSE 9,115 Mapelli and Giacobbo (2018) 
Isolated binary population synthesis, MOBSE 78 Artale et al. (2019) 
Isolated binary population synthesis, MOBSE 4, 37] Baibhav et al. (2019) 
Isolated binary population synthesis, SEVN 20 Boco et al. (2019) 
Isolated binary population synthesis, StarTrack 5, 230 Chruslinska et al. (2019) 
Isolated binary population synthesis, BPASS 209, 269] Eldridge et al. (2019) 
Isolated binary population synthesis, COMPAS 19, 204] Neijssel et al. (2019) 
Isolated binary population synthesis, StarTrack 0.48, 297] Belczynski et al. (2020) 
Isolated binary population synthesis, MOBSE 6, 80] Giacobbo and Mapelli (2020) 
Isolated binary population synthesis, MOBSE 49*35 Santoliquido et al. (2020) 
Isolated binary population synthesis, BPASS 58, 6225] Tang et al. (2020) 
Isolated binary population synthesis, COSMIC 3.7, 1100] Zevin et al. (2020) 
Isolated binary population synthesis, BPASS 8.7, 498] Ghodla et al. (2021) 
Isolated binary population synthesis, StarTrack 4, 16] Olejak et al. (2021) 
Isolated binary population synthesis, POSYDON 5.7, 77] Román-Garza et al. (2021) 
Isolated binary population synthesis, MOBSE 1.8,128] Santoliquido et al. (2021) 
Isolated binary population synthesis, BSE 10, 72] Shao and Li (2021) 
Isolated binary population synthesis, BPASS, 27 Briel et al. (2022b) 
Isolated binary population synthesis, COMPAS 2.2, 830] Broekgaarden et al. (2021, 2022) 
Isolated binary population synthesis, StarTrack 4, 27] Olejak et al. (2022) 
Chemically homogeneous evolution, MESA 0.02, 0.2] Marchant et al. (2017) 
Population III stars 0.0002, 0.016] Belczynski et al. (2017) 
Hierarchical triples 1.9 x 10-7, 22]  Fragione and Loeb (2019a,b) 
Hierarchical triples, SecularMultiple 345, 680] Hamers and Thompson (2019) 
Hierarchical triples in young star clusters, OKINAMI, 0.04, 0.34 Trani et al. (2021) 
Hierarchical quadruples, MSE 5.7, 57.1] Vynatheya and Hamers (2021) 
Wide isolated binaries with flybys 40 Michaely and Naoz (2022) 
Globular cluster dynamics 0.01, 0.12 Clausen et al. (2013) 
Globular cluster dynamics, ARCHAIN $0.1 Arca Sedda (2020b) 
Globular cluster dynamics, CMC 0.009, < 5.5] Ye et al. (2020) 
Nuclear star cluster dynamics, with SMBH 0.02, 0.4] Petrovich and Antonini (2017) 
Nuclear star cluster dynamics, with SMBH 0.17, 0.52 Fragione et al. (2019) 
Nuclear star cluster dynamics, with SMBH = [2,5] Stephan et al. (2019) 
Nuclear star cluster dynamics, S 0.01 Arca Sedda (2020b) 
Nuclear star cluster dynamics, with SMBH < 300 McKernan et al. (2020) 
Nuclear star cluster dynamics, with SMBH = [0.15, 0.3] Wang et al. (2021) 
Young star clusters, starlab S04 Ziosi et al. (2014) 
Young/Open star clusters, Nbody7 <3 x 10-3 Fragione and Banerjee (2020) 
Young star clusters, MOBSE S 28 Rastello et al. (2020) 
Young star clusters, MOBSE 4135 Santoliquido et al. (2020) 
Young star clusters (0.04, 36.6] Arca Sedda (2021) 


Rates of compact object coalescences 9 


3.1 Gravitational waves 


Gravitational-wave data should, at first glance, enable straightforward rate 
estimates. These are direct observations of compact-object coalescences. The 
data are homogeneous, with well-understood selection effects which can be 
modelled by injecting mock signals into the data stream and measuring how 
many of them will be successfully extracted by the search pipelines. Once the 
surveyed timex volume of the instruments is thus estimated, the volumetric 
merger rate is given by the number of detections in this timex volume. We 
might therefore expect that analysing some 70 confident BH-BH mergers from 
the first, second, and third advanced detector observing run in Abbott et al. 
(2021e) would enable the BH-BH merger rate to be inferred to a fractional 
accuracy of 1/4/70 ~ 12% based on Poisson statistics. 

However, there are several complicating factors to this story, as a result of 
which the total binary black hole merger rate is still uncertain at the factor of 
two level. Firstly, the detection rate is a sensitive function of mass, with more 
massive binaries typically providing higher-amplitude gravitational-wave sig- 
nals that can be detected through a greater volume. Therefore, inferring the 
total merger rate relies on simultaneously measuring the mass distribution of 
merging black holes, which requires many more observations (Abbott et al., 
2021e). One way to partially overcome this uncertainty is to focus on observa- 
tions in a more narrow mass band (Roulet et al., 2020; Abbott et al., 2021c); 
for example, Roulet et al. (2020) find that the merger rate of 20-30 Mo BHs? 
is 1.5-5.3 Gpc ? yr-!, while Abbott et al. (2021c) find a BH-BH merger rate 
in the same mass range of 3.7* 1$ or 2.6* 12 Gpc ? yr-! depending on the as- 
sumed mass distributions. The inferred merger rates reported in the summary 
tables and figures as “PDB (ind)”, “MS”, and “BGP” refer to the correspond- 
ing mass distribution models described in Abbott et al. (2021e): Power law + 
dip 4- break (with independent mass pairings for the two components), Power 
law + peak multi-source (separate BH mass models in BH-BH and NS-BH 
binaries) and Binned Gaussian process (unmodelled with regularisation). 

Secondly, another uncertainty arises from the possibility that the merger 
rate varies with redshift, for which there is already significant support in the 
data (Fishbach et al., 2021; Abbott et al., 2021e). Allowing for the variation of 
the merger rate with redshift (but not yet of the mass distribution with red- 
shift, although some theoretical models predict this) further shifts the binary 
black hole merger rate at redshift z = 0. The “z-dependent” rate reported for 
BH-BH mergers comes from the redshift-dependent merger rate model evalu- 
ated at z — 0 (Abbott et al., 2021e). 

Thirdly, details of search methodology and tuning can impact the set of 
detections (e.g., Abbott et al., 2021b,d), and independent searches over the 
publicly released LIGO-Virgo data have yielded slightly different source popu- 
lations (e.g., Zackay et al., 2019; Venumadhav et al., 2020; Nitz et al., 2021a,b). 
While searches with slightly different sensitivities and different numbers of de- 


2 We generally use solar masses (Mc) and solar radii (Rc) as mass and distance units. 
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Fig. 2: NS-BH merger rates from Table 2. Notation as in Figure 1. (GitHub) 


tections should still yield consistent rate estimates, some fluctuation is to be 
expected, especially through the impact of different inferred mass distributions 
as described above. 

Fourthly, while the chirp mass m2/?m2/? (m, + ma)- !/? is measured accu- 
rately for low-mass binaries through gravitational-wave signals, the mass ratio 
measurement is much less accurate, making it difficult to distinguish neutron 
stars from low-mass black holes (see, e.g., Hannam et al. 2013). A measurement 
of tidal deformation could point unambiguously to a neutron star rather than 
a black hole, but such measurements require sensitivity at higher frequencies, 
and only upper limits have been placed so far (Abbott et al., 2017b). There- 
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fore, in the absence of an electromagnetic counterpart, it is not clear how to 
classify a system such as GW190814 (Abbott et al., 2020d), which consisted 
of a ~ 23 Mo BH coalescing with a ~ 2.6 Me compact object that is probably 
a low-mass black hole but could be a very high-mass neutron star. 

Estimates of the NS-NS and NS-BH merger rates are based on a very small 
number of confident detections, including GW170817 and GW190425 in the 
former category and GW200105 and GW200115 in the latter category (Abbott 
et al, 2017b, 2020a, 2021a). Furthermore, these estimates are particularly 
sensitive to assumptions about the poorly constrained mass distributions in 
NS-NS and NS-BH binaries. 


3.2 Short gamma-ray bursts and kilonovae 


Mergers of two neutron stars have long been expected to be accompanied 
by a short burst of spectrally hard gamma radiation accompanying an ultra- 
relativistic jet, with subsequent synchrotron radiation afterglow from the in- 
teraction of the jet with the circumstellar medium (Paczyński and Rhoads, 
1993; Rosswog and Ramirez-Ruiz, 2002; Nakar, 2007). Meanwhile, the for- 
mation and subsequent radioactive decay of neutron-rich elements through 
r-process nucleosynthesis was conjectured to power an optical kilonova (Li 
and Paczyński, 1998; Rosswog et al., 1999; Metzger et al., 2010; Kasen et al., 
2013). Both expectations were spectacularly confirmed with the detection of 
a short gamma-ray burst (Abbott et al., 2017a), followed by the observation 
of a kilonova (Abbott et al., 2017c; Kasen et al., 2017; Metzger, 2019), and 
finally a broad-spectrum (radio through optical to X-ray) afterglow (Mooley 
et al., 2018; Lyman et al., 2018; Troja et al., 2017) following the NS-NS merger 
GW17017. 

One concern with using short gamma-ray bursts (SGRBs) to infer binary 
neutron star coalescence rates is the uncertainty in a one-to-one association 
between these events. On the one hand, some SGRBs and perhaps even kilono- 
vae could be associated with NS-BH mergers (Troja et al., 2008; Berger, 2014; 
Li et al., 2017; Gompertz et al., 2020). However, this is likely only possible for 
those events when the neutron star is tidally disrupted at a sufficient distance 
from the black hole, requiring the black hole to be low in mass and/or rapidly 
rotating (Foucart et al., 2018). Therefore, the contribution of NS-BH mergers 
to SGRBs is likely low (e.g., Drozda et al., 2020). On the other hand, given the 
challenges in reproducing ultra-relativistic jet formation in numerical merger 
simulations (Mósta et al., 2020), it cannot be determined with confidence 
whether all neutron star mergers should be accompanied by SGRBs. 

Nevertheless, SGRBs have been used for more than a decade as proxies for 
the NS-NS merger rate, with Nakar et al. (2006) estimating the merger rate as 
at least 10 Gpc^? yr^ 1, but possibly several orders of magnitude higher. The 
key challenge is the uncertainty in the angular size of the jet: the smaller the 
jet, the smaller fraction of SGRBs beamed toward Earth, and hence the larger 
the intrinsic rate based on the observed rate of SGRBs. Fong et al. (2012) 
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Table 3: Summary of local BH-BH coalescence rates, part 1: observations and 
isolated binary evolution model predictions (Supplementary material and Zen- 


odo) 
Source Rate Reference 

[Gpc7? yr] 
Observations: 
Gravitational-wave observations, PDB (ind) 2277 GWTC-3; Abbott et al. (2021e) 
Gravitational-wave observations, MS arcs GW'TC-3; Abbott et al. (2021e) 
Gravitational-wave observations, BGP 850718 0 GW'TC-3; Abbott et al. (2021e) 
Gravitational-wave observations, z-dependent artus. GW'TC-3; Abbott et al. (2021e) 
Models: Rate 
Isolated binary population synthesis, Star Track 2, 40] O'Shaughnessy et al. (2010) 
Isolated binary population synthesis, Brussels code « 96, 1140] Mennekens and Vanbeveren (2014) 
Isolated binary population synthesis, Star Track 14, 2500] de Mink and Belczynski (2015) 
Isolated binary population synthesis, Star'Track 0.5, 221] Dominik et al. (2015) 
Isolated binary population synthesis, BSE 850 Lamberts et al. (2016) 
Isolated binary population synthesis, Scenario Machine 100 Lipunov et al. (2017) 
Isolated binary population synthesis, MOBSE 20,572] Mapelli et al. (2017) 
Isolated binary population synthesis, BSE 20, 320] Ablimit and Maeda (2018) 
Isolated binary population synthesis, Star Track 32, 1072 Chruslinska et al. (2018) 
Isolated binary population synthesis, MOBSE 43, 1500 Giacobbo and Mapelli (2018) 
Isolated binary population synthesis, Star' Track 89, 203] Klencki et al. (2018) 
Isolated binary population synthesis, COMBINE 0.6, 109] Kruckow et al. (2018) 
Isolated binary population synthesis, MOBSE 146, 240 Mapelli and Giacobbo (2018) 
Isolated binary population synthesis, MOBSE 142 Artale et al. (2019) 
Isolated binary population synthesis, MOBSE 30, 60] Baibhav et al. (2019) 
Isolated binary population synthesis, Star Track 12,1072 Chruslinska et al. (2019) 
Isolated binary population synthesis, BPASS 56, 174]. Eldridge et al. (2019) 
Isolated binary population synthesis, COMPAS 59, 1157 Neijssel et al. (2019) 
Isolated binary population synthesis, SEVN 90 Spera et al. (2019) 
Isolated binary population synthesis, StarTrack 1.24, 1368] Belczynski et al. (2020) 
Isolated binary population synthesis, MOBSE 43, 160 Giacobbo and Mapelli (2020) 
Isolated binary population synthesis, MOBSE 6, 37] Mapelli et al. (2021) 
Isolated binary population synthesis, COMPAS 51,87] Riley et al. (2021) 
Isolated binary population synthesis, POSYDON 70, 203 Román-Garza et al. (2021) 
Isolated binary population synthesis, MOBSE 50732 Santoliquido et al. (2020) 
Isolated binary population synthesis, BPASS 10, 219 Tang et al. (2020) 
Isolated binary population synthesis, COSMIC 84, 6900] Zevin et al. (2020) 
Isolated binary population synthesis, POSYDON 39, 170 Bavera et al. (2021) 
Isolated binary population synthesis, BPASS 31,873 Ghodla et al. (2021) 
Isolated binary population synthesis, Star' Track 18, 89] Olejak et al. (2021) 
Isolated binary population synthesis, MOBSE 10, 105.4] Santoliquido et al. (2021) 
Isolated binary population synthesis, BSE 43, 76] Shao and Li (2021) 
Isolated binary population synthesis, BPASS, 6.5 Briel et al. (2022b) 
Isolated binary population synthesis, COMPAS 3.8,810] Broekgaarden et al. (2021, 2022) 
Isolated binary population synthesis, SeBa 19.1, 101.7] Dorozsmai and Toonen (2022) 
Isolated binary population synthesis, Star' Track 43, 102] Olejak et al. (2022) 
Chemically homogeneous evolution 2, 80] Mandel and de Mink (2016) 
Chemically homogeneous evolution, MESA 0.7, 16] Marchant et al. (2016) 
Chemically homogeneous evolution, MESA 5.8, 7] du Buisson et al. (2020) 
Chemically homogeneous evolution, COMPAS 4, 32] Riley et al. (2021) 
Population III stars 12, 25] Kinugawa et al. (2014) 
Population III stars, StarTrack 0.016, 1.9] Belczynski et al. (2017) 
Population III stars, BSE 0.38, 2.9] Hijikawa et al. (2021) 
Population III stars, BSE 3.34, 21.2] Kinugawa et al. (2020) 
Population III stars in Nuclear Star Clusters 0.02, 0.6] Liu and Bromm (2021) 
Population III stars, BSE 0.1 Tanikawa et al. (2021) 
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Table 3: Summary of local BH-BH coalescence rates, part 2: dynamical for- 
mation models and primordial black hole predictions 


Source Rate Reference 

Gpc^? yr-1] 
Hierarchical triples, Rebound package 0.14, 6.3] Silsbee and Tremaine (2017) 
Hierarchical triples, TRES 0.3, 2.5] Antonini et al. (2017) 
Hierarchical triples 2, 25] Rodriguez and Antonini (2018) 
Hierarchical triples in globular clusters, CMC = 0.35, 1] Martinez et al. (2020) 
Hierarchical triples in young star clusters 0.2, 1.44] Trani et al. (2021) 
Hierarchical quadruples, MSE 8.5, 39.8] Vynatheya and Hamers (2021) 
Hierarchical triples, MOBSE and TSE 12.5, 95.2] Stegmann et al. (2022) 
Wide isolated binaries with flybys, COMPAS 1,20 Raveh et al. (2022) 
Globular cluster dynamics 7.25, 29] Bae et al. (2014) 
Globular cluster dynamics, CMC 3.8,13 Rodriguez et al. (2015) 
Globular cluster dynamics 5 Antonini and Rasio (2016) 
Globular cluster dynamics, CMC 2,20 Rodriguez et al. (2016a) 
Globular cluster dynamics, MOCCA 5.4, 30 Askar et al. (2017) 
Globular cluster dynamics 13, 57] Fujii et al. (2017) 
Globular cluster dynamics, Nbody6 6.5, 26 Park et al. (2017) 
Globular cluster dynamics 4, 60 Fragione and Kocsis (2018) 
Globular cluster dynamics 0.8, 20 Hong et al. (2018) 
Globular cluster dynamics, CMC 4,18 Rodriguez and Loeb (2018) 
Globular cluster dynamics 26 Choksi et al. (2019) 
Globular cluster dynamics, cBHBd 0.2, 50 Antonini and Gieles (2020) 
Globular cluster dynamics, CMC 9, 30 Kremer et al. (2020) 
Globular cluster dynamics, MOBSE 0.8, 7] Mapelli et al. (2021) 
Nuclear star cluster dynamics, without SMBH 1,10 Miller and Lauburg (2009) 
Nuclear star cluster dynamics, with SMBH 0.002, 0.6] Antonini and Perets (2012) 
Nuclear star cluster dynamics, without SMBH 1.5 Antonini and Rasio (2016) 
Nuclear star cluster dynamics:, with SMBH ~ 1.2 Bartos et al. (2017) 
Nuclear star cluster dynamics, with SMBH 0.6, 15] Petrovich and Antonini (2017) 
Nuclear star cluster dynamics, with SMBH ~3 Stone et al. (2017) 
Nuclear star cluster dynamics, with SMBH 0.01, 0.4] Hamers et al. (2018) 
Nuclear star cluster dynamics, with SMBH 1, 3] Hoang et al. (2018) 
Nuclear star cluster dynamics, with SMBH 0.17, 0.52] Fragione et al. (2019) 
Nuclear star cluster dynamics, with SMBH 0.002, 0.04] Rasskazov and Kocsis (2019) 
Nuclear star cluster dynamics, with SMBH = [7,15] Stephan et al. (2019) 
Nuclear star cluster dynamics with SMBH 3.3, 8.6] Arca Sedda (2020a) 
Nuclear star cluster dynamics with SMBH 0.002, 18] Gróbner et al. (2020) 
Nuclear star cluster dynamics, with SMBH 0.02, 60] Tagawa et al. (2020) 
Nuclear star cluster dynamics, with SMBH 0.1, 1.6]. Yang et al. (2020) 
Nuclear star cluster dynamics, with SMBH 6, 20] Ford and McKernan (2021) 
Nuclear star cluster dynamics 0.09, 10] Mapelli et al. (2021) 
Nuclear star cluster dynamics, with SMBH = [0.3, 5] Wang et al. (2021) 
Nuclear star cluster dynamics, with SMBH 0.083, 14] Boekholt et al. (2022) 
Nuclear star cluster dynamics S1 Codazzo et al. (2022) 
Young star clusters, starlab S15 Ziosi et al. (2014) 
Young star clusters, starlab 21 Mapelli (2016) 
Open star clusters, NBODY7 <2 Rastello et al. (2019) 
Young star clusters, NBODY6 < [55, 110] Di Carlo et al. (2020) 
Open star clusters, NBODY6 [35, 70] Kumamoto et al. (2020) 
Young star clusters, MOBSE pium Santoliquido et al. (2020) 
Young and open star clusters, NBODY7 [0.5, 37.9] Banerjee (2021) 
Young star clusters, MOBSE [0.1, 18] Mapelli et al. (2021) 
Young star clusters, NBODY6 MOBSE BS. Rastello et al. (2021) 
Primordial binaries [0.02, < 3] Bird et al. (2016) 
Primordial binaries < [0.2, 105] Ali-Haimoud et al. (2017) 
Primordial binaries < [6, 104] Raidal et al. (2019) 
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Fig. 3: BH-BH merger rates from Table 3. Notation as in Figure 1. (GitHub) 
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used evidence for jet breaks to constrain the angular size of the jet, allowing 
them to infer an NS-NS merger rate. However, the evidence from GW170817 
firmly indicated that the jet has spatial structure, rather than a simple top- 
hat profile (Troja et al., 2018; Resmi et al., 2018; Lamb et al., 2019). This 
requires a re-evaluation of earlier models, an ongoing effort (Paul, 2018); for 
example, Margutti and Chornock (2020, see their figure 3) estimate the rate 
of GW170817A-like SGRBs as > 300 Gpc ? yr-! given the observed off-axis 
angle under the assumption that other SGRBs would have similar jet profiles. 

The first kilonova observation was made in 2013 (Tanvir et al., 2013). 
While more candidates have been observed in the last few years, they are 
still difficult to unambiguously distinguish from other transients. At present, 
their luminosity function is not yet sufficiently well known to infer an accurate 
rate. However, the absence of transients with a kilonova signal similar to that 
accompanying GW170187 in data from the Palomar Transient Factor allowed 
Kasliwal et al. (2017) to place a conservative 3-0 upper limit of 800 Gpc^? 
yr-!. Andreoni et al. (2021, see their Fig. 9) compile a list of optical kilonova 


surveys with upper limits between 900 and 99000 Gpc ? yr7!. 


3.3 Pulsars in Galactic compact object binaries 


Observations of double compact object binaries prior to merger are the next 
best thing to direct merger observations for measuring merger rates. Pro- 
vided that the masses, period or separation, and eccentricity of the binary are 
known, the time to merger through the emission of gravitational waves can 
be computed (Peters, 1964). Assuming a steady-state configuration, which is 
a reasonable assumption for our Galaxy with its slowly evolving metallicity 
and star formation rate, an observed population of pre-merger compact object 
binaries yields an estimate of the coalescence rate. 

So far, NS-NS binaries are the only double compact objects for which direct 
observations are available. Nearly 20 such binaries have been observed in the 
Galaxy through radio pulsar pulsar observations (see, e.g., Tauris et al., 2017; 
Farrow et al., 2019, for recent compilations). Of course, these are expected to 
represent a relatively small fraction of all Galactic double neutron stars, with 
most remaining unobserved because neither neutron star is a pulsar or any 
pulsars are too dim or beamed away from the Earth. Moreover, pulsars in the 
tightest binaries are rare because of their rapid inspiral and may be challenging 
to detect because the pulsations are Doppler shifted by orbital modulation. 
'Therefore, the key challenge in inferring the merger rate of close double neutron 
stars in the Galaxy from these observations lies in accounting for selection 
effects, particularly the uncertain luminosity function and beaming fraction of 
pulsars. 

Fortunately, the rapidly growing amount of pulsar data since the pioneering 
efforts of Phinney (1991) and Narayan et al. (1991) to infer NS-NS merger 
rates from binary pulsar observations enables increasingly accurate treatments 
of these selection effects. Much of the recent literature follows the statistical 
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framework of Kim et al. (2003), with the addition of varying models of the 
luminosity distribution (Kalogera et al., 2004), more sophisticated models of 
the pulsar beaming fraction which could vary with pulsar age (O'Shaughnessy 
and Kim, 2010), an attempt to account for the observation of both pulsars 
in the double pulsar J0737-3039 (Kim et al., 2015) and the addition of new 
systems (with the relatively small number of known Galactic binaries, a single 
new short-lived double neutron star could impact the rate estimate at the 
factor of two level, as shown by Kim et al. 2006). We report the rates from a 
few recent papers to indicate the evolving estimates of the double neutron star 
merger rate (see Abadie et al., 2010, for a compilation including older results). 

'There are several complicating factors in extracting coalescence rates from 
the observed Galactic double neutron stars. B2127+11C and possibly J1807- 
2500 are known double neutron stars in Galactic globular clusters; these are 
generally excluded from rate estimates because they are believed to be rare 
dynamically formed systems (Ye et al., 2020, but see Andrews and Mandel 
2019). Another concern is that the detectability of a Galactic double neu- 
tron star may correlate with the evolutionary channel and with the system 
properties, e.g., through the amount of recycling the pulsar experiences. This 
could lead some sub-populations of pulsars — perhaps the more massive ones 
like GW190425 (Abbott et al., 2020a; Romero-Shaw et al., 2020; Vigna-Gómez 
et al., 2021a) — to be missed in Galactic surveys. Moreover, extrapolation from 
the Milky Way merging double neutron star rate to a volumetric rate depends 
on the accuracy of estimating the local star formation rate (de Freitas Pacheco 
et al., 2006) and on the sensitivity of the yield of merging double neutron stars 
to metallicity. 

No Galactic radio pulsars in NS-BH binaries have been detected so far. 
However, the lack of such observations only weakly constrains the NS-BH 
merger rate from above. If the heavier BH forms first, which appears intuitive 
given the shorter lifetimes of more massive stars, the NS will not be recycled 
by accretion, which limits the time during which it can be observed as a radio 
pulsar (Chattopadhyay et al., 2021). Mass transfer can invert the mass ratio, 
allowing the NS to form first from the originally more massive star; however, 
the formation rate of such systems is likely quite low (Pfahl et al., 2005). 


3.4 High-mass X-ray binaries 


As we will discuss in the next section, most compact object binaries formed 
through isolated binary evolution are expected to go through a high-mass X- 
ray binary stage, in which the first companion to form a black hole or neutron 
star accretes from the wind of the stellar companion. This phase may directly 
precede the formation of two compact objects, and can therefore be used to 
constrain compact binary coalescence rates. For example, Bulik et al. (2011) 
used the extragalactic high-mass BH X-ray binaries IC10 X-1 and NGC300 
X-1 to estimate the merger rate of BH-BH binaries at 360200 Gpc ? yr“. 
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However, such estimates suffer from significant uncertainties. Firstly, there 
are often large systematic errors in the masses and separations of the observed 
systems; in particular, the masses of IC10 X-1 and NGC300 X-1 may be much 
lower than in the analysis of Bulik et al. (2011) because orbital velocities 
are very challenging to measure in the presence of optically thick Wolf-Rayet 
winds (Laycock et al., 2015). Even well-studied high-mass X-ray binaries such 
as Cygnus X-1 (Orosz et al., 2011) may require significant corrections to their 
masses (Miller-Jones et al., 2021), while the uncertain present properties of 
Cygnus X-3 translate into uncertainty about its future fate (Belczynski et al., 
2013). Secondly, the observed systems typically come from a heterogeneous set 
of observations, complicating the determination of selection effects. Thirdly, 
the future evolution of such systems is sensitive to assumptions about the 
dynamical stability of any future mass transfer episodes (van den Heuvel et al., 
2017), the amount of mass loss in winds (cf. Belczynski et al. 2011 and Neijssel 
et al. 2021 for the case of Cygnus X-1), and the remnant mass and natal kick 
of the remaining supernova. 

Inoue et al. (2016) and Finke and Razzaque (2017) use observations of 
ultra-luminous X-ray sources (ULXs) with X-ray luminosities exceeding 10?? erg s7 
to estimate BH-BH merger rates of € 100 Gpc ? yr-! and 15-400 Gpc ? yr-!, 
respectively. These estimates, however, are highly uncertain. They rely on the 
assumption that all ULXs are powered by accretion onto a BH, which is contra- 
dicted by observations that at least some ULXs (perhaps all with confidently 
identified compact objects?) are powered by accretion onto an NS (Bachetti 
et al., 2014). They also make simplifying assumptions about the future evolu- 
tionary history of ULXs. 


1 


3.5 Other observational constraints 


Abadie et al. (2010) suggested that a strict upper limit of < 5 x 10* mergers 
Gpc^? yr-! could be placed on the NS-NS merger rate by assuming that at 
least one companion must have been formed through a type Ib/Ic supernova 
(taking an optimistic rate estimate from Cappellaro et al. 1999). However, this 
upper limit is likely a very significant over-estimate, as most type Ib/Ic super- 
novae happen in other systems than merging binary neutron stars (perhaps 
even in systems that had already been disrupted, Hirai et al. 2020). At the same 
time, the type Ib/Ic supernova rate cannot even be considered a strict upper 
limit, because many merging neutron star binaries could form through a com- 
bination of electron-capture supernovae and ultra-stripped supernovae rather 
than classical core-collapse supernovae (Podsiadlowski et al., 2004; Dall'Osso 
et al., 2014; Tauris et al., 2015, 2017; Vigna-Gómez et al., 2018), and some 
compact remnants may form in the absence of any explosions. 

Signatures of chemical enrichment associated with binary neutron star 
mergers, particularly the evolution of r-process nucleosynthesis elements, could 
serve as a useful constraint on the binary neutron star merger rate (Men- 
nekens and Vanbeveren, 2014; Hotokezaka et al., 2015; Beniamini et al., 2016). 
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However, Siegel et al. (2019); Kobayashi et al. (2020) propose that magneto- 
rotational supernovae could instead be responsible for r-process enrichment. 

Compact-object mergers or the subsequent evolution of compact-object 
merger products have been proposed as a source for some fast radio bursts 
(Totani, 2013; Ravi and Lasky, 2014); however, there are likely multiple pro- 
genitor channels for fast radio bursts (Pleunis et al., 2021), so these cannot 
presently be confidently used for merger rate estimates. 


4 Modelled rates 


In this section, we summarise some of the ongoing efforts to model the rates 
of compact object mergers formed through a variety of evolutionary channels. 
'These channels offer alternative pathways to overcoming a key challenge in 
gravitational-wave astronomy. Gravitational-wave emission is only efficient in 
close binaries: in order for a circular binary with equal-mass components of 
mass M to merge within a time T' exclusively through radiation reaction from 
gravitational-wave emission, the separation a must be smaller than 


M \3/4 T 1/4 
<3. = — ; 1 
“Sart (3) (wom) o 


Thus, to merge within the current age of the Universe, 14 Gyr, two 1.4 Mo 
neutron stars must be fewer than 5 solar radii apart, while even two 35 Mo 
black holes must be within a quarter of an astronomical unit. 

On the other hand, massive stars typically expand to several au (hundreds 
to thousands of solar radii; 1 au ~ 215 Ro) in size during their evolution. 
Therefore, the key challenge to forming a merging compact-object binary lies 
in fitting a peg into a hole that is one or two orders of magnitude smaller than 
the size of the peg. The proposed solutions generally take the form of either 
relying on stellar evolution and mass transfer to overcome the problem in an 
isolated binary, or circumventing the problem altogether by relying on dy- 
namics to form binaries with a sufficiently small separation that gravitational- 
wave emission can drive them to merge. Here we sketch out the key steps in 
the proposed formation channels, undertake back-of-the-envelope estimates of 
the rates, following Mandel and Farmer (2022), and describe some of the key 
uncertainties leading to the typically broad range of model predictions. 

One particular source of uncertainty that is common to all models described 
below is the uncertain rate and metallicity of star formation across cosmic his- 
tory. Even binaries merging in the relatively local Universe (detections through 
the third LIGO/Virgo observing run had moderate redshifts z < 1, Abbott 
et al. 2021d) could have formed at high redshift (Belczynski et al., 2016). 
Therefore, the rate of star formation at higher redshift influences local merger 
rates. Moreover, the yield of coalescing compact objects per unit star forming 
mass depends strongly on metallicity (Dominik et al., 2013; Neijssel et al., 
2019). This is particularly true for black holes, whose massive progenitors can 
lose large fractions of their mass in winds in high-metallicity environments, 
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widening the binaries and reducing merger rates. Therefore, the local merger 
rate is sensitive to the poorly known metallicity-specific star formation rate 
at high redshifts (Belczynski et al., 2010; Mapelli et al., 2017; Klencki et al., 
2018; Chruslinska et al., 2019; Neijssel et al., 2019; Chruslinska and Nelemans, 
2019; Tang et al., 2020; Briel et al., 2022a; Broekgaarden et al., 2021, 2022; 
Santoliquido et al., 2021; van Son et al., 2022). 


4.1 Isolated binary evolution modelling (population synthesis) 


Almost all stars massive enough to form neutron stars and black holes are born 
in binaries or higher-multiplicity systems (Sana et al., 2012; Moe and Di Ste- 
fano, 2017). We begin by considering compact-object coalescences in isolated 
binaries, namely those which do not appreciably interact with their environ- 
ment. This channel has been studied since the 1970s. Some of the ground- 
breaking works particularly relevant to the formation of coalescing compact- 
object binaries include Tutukov and Yungelson (1973); van den Heuvel and De 
Loore (1973); De Loore et al. (1975); Flannery and van den Heuvel (1975); Tu- 
tukov and Yungelson (1993); Lipunov et al. (1997); Bethe and Brown (1998); 
Dewi et al. (2006); Nelemans et al. (2001); Voss and Tauris (2003); Kalogera 
et al. (2007); O’Shaughnessy et al. (2008) (see Mandel and Farmer 2022 for a 
brief summary of the early history). 


4.1.1 Isolated binary evolution with mass transfer 


In this channel, two massive stars are born in a wide binary, giving them 
sufficient space to evolve. The binary shrinks as a consequence of mass transfer 
just when the stars themselves shrink in radius, allowing gravitational waves 
to take over. 

We begin by sketching out a very simplified version of binary evolution 
in this model. The primary - the initially more massive star — completes core 
hydrogen fusion first. At the end of the main sequence, its helium-rich core con- 
tracts and the hydrogen-rich envelope expands. As the star fills the so-called 
Roche lobe, tidal gravity from the companion overcomes the self-gravity of 
the primary, and mass transfer onto the secondary commences. Eventually, 
the primary leaves behind a naked helium star, which continues nuclear fu- 
sion, losing mass through winds, until it collapses into a compact object. The 
mass loss and/or the natal kick from asymmetric mass ejection accompanying 
the supernova may unbind the binary. If it does not, the binary continues its 
evolution until the secondary expands and commences mass transfer onto the 
now compact-object primary. Because of the mass gain of the secondary during 
the first episode of mass transfer and the mass loss of the primary due to mass 
transfer, winds, and the supernova, the secondary may be significantly more 
massive than the primary at this stage. Consequently, mass transfer could 
harden the binary (van den Heuvel et al., 2017), and could even become dy- 
namically unstable, leading to the formation of a common envelope (Paczyński, 
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1976). Once the common envelope is not co-rotating with the binary, drag 
forces dissipate orbital energy, which may allow for common-envelope ejection 
after the binary hardens by two or three orders of magnitude (Ivanova et al., 
2013). At this stage, the secondary has lost its envelope to become a naked 
helium star, so it can fit into the tight binary. (It is also possible that both 
stars are evolved giants during the first mass transfer interaction, leading to 
the simultaneous ejection of both envelopes during a single common-envelope 
event, as described by Dewi et al. 2006.) After more mass loss through winds 
the secondary may also collapse into a compact object. Assuming the second 
supernova leaves behind a bound, compact binary, this binary will eventually 
merge by emitting gravitational waves. 

We can crudely estimate the expected merger rate through this channel 
with a Drake-like equation describing the fraction of all binaries that go on 
to make merging double compact objects of a particular type (Mandel and 
Farmer, 2022), with chosen values corresponding to binary black hole forma- 
tion: 


fpco = fprimary x fsecondary x finit sep x feurvive SN1 X fcr x dourvive SN2 
X fmerge ~ 0.001 x 0.5 x 0.5x 1x 0.1x 1x 0.225 x 107°. (2) 


Assuming a local star formation rate of ~ 10" Mz Gpc^? yr-! (Madau and 
Dickinson, 2014) with an average mass of ~ Mo, a yield of fpco = 5 x 1079 
would lead to a compact-object merger rate of ~ 50 Gpc^? yr7!. 

Here we describe the terms in this equation, which allows us to highlight the 
broad range of physics involved and to comment on the key uncertainties. In 
practice, most modern models rely on some form of population synthesis, where 
large numbers of stellar binaries are evolved in order to compute the yield, 
frequently by using simplified recipes to reduce computational cost (Postnov 
and Yungelson, 2014; Han et al., 2020). 

Stars with initial mass at zero age on the main sequence of roughly > 20 Mo 
can form a BH, so the fraction of all stars drawn from the initial mass function 
that will have mass in this range is fprimary © 0.001. Meanwhile, for forming an 
NS, the zero-age main sequence mass should be roughly between 8 and 20 Mo, 
SO fprimary % 0.002. The mass ratio between secondary and primary is roughly 
uniformly distributed (Sana et al., 2012), so if the primary falls in the mass 
range of interest, roughly half the secondaries will do so, i.e., fsecondary © 0.5. 
The fraction of binaries with initial separations in the range to allow for the 
first episode of mass transfer from evolved donors is quite large, finit sep © 0.5, 
given the significant amount of expansion during this phase and the uniform 
in the logarithm distribution of initial separations (Opik, 1924). 

The probabilities of the binary surviving the first and second supernovae 
without being disrupted, feurvive sN1 and feurvive SN2, are both close to 1 for the 
heavier black holes that may form through complete collapse with negligible 
natal kicks (Mirabel, 2017, but see Repetto et al. 2012; Mandel 2016; Atri et al. 
2019). However, the high typical natal kicks of neutron stars, of order 300 km 
s^! (Hobbs et al., 2005; Verbunt et al., 2017), can disrupt > 95% of wide bina- 
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ries in the first supernova (Vigna-Gómez et al., 2018; Renzo et al., 2019). Re- 
duced kicks in electron-capture supernova (Podsiadlowski et al., 2004; Gessner 
and Janka, 2018) may protect some of the wide binaries from unbinding, in- 
creasing fsurvive SN1- By the time of the second supernova, the binary will have 
gone through a common envelope in this channel, hardening it and making 
disruption less likely; moreover, low-mass helium star progenitors of neutron 
stars may re-expand after the helium main sequence, engage in another mass 
transfer episode (case BB mass transfer) and ultimately explode in very weak 
ultra-stripped supernovae with greatly reduced kicks (Tauris et al., 2015, but 
see Mandel and Müller 2020). Therefore, it is plausible that fiurvive sno 7 1, 
even for double neutron stars, although observed eccentricities of ~ 0.6 for 
a subset of Galactic double neutron stars, including the Hulse- Taylor binary, 
could point to a subpopulation of systems that avoid ultra-stripping and may 
be disrupted by the second supernova. 

'The common-envelope phase is perhaps the least certain one in massive 
binary evolution (Dominik et al., 2012). There are questions about both the 
onset of the common-envelope phase (do the mass ratio and the evolutionary 
state of the donor and accretor make mass transfer dynamically unstable?) and 
about the ultimate outcome of the evolution (is sufficient energy deposited in 
the envelope to eject it, or does this episode end in merger?). Typical popu- 
lation synthesis models based on current, very simplified recipes for both the 
stability mass transfer and the survivability of the common-envelope phase 
yield for 0.1. 

The fraction of compact-object binaries that ultimately coalesce in the age 
of the Universe fmerge is set by the distribution of masses and separations. If 
the post-common-envelope binary separation is distributed uniformly in the 
logarithm, the delay time r between formation and merger, which scales with 
the fourth power of the separation, will also follow a log-uniform distribution, 
i.e., p(T) x 1/7. In that case, the logarithmic distribution allows for typical 
fmerge ~ 0.2; however, as mentioned above, this can be significantly reduced, 
e.g., if winds increase the BH binary separation at high metallicity. 

We now describe the key uncertainties in these estimates. Different treat- 
ments of some of the physics summarised below are responsible for the broad 
range of predictions reported in the executive summary. 

'The first significant uncertainty lies in the distribution of initial conditions: 
the masses of the two companions, their separations, and the binary eccentric- 
ity (Sana et al., 2012; Moe and Di Stefano, 2017). These alone can impact the 
merger rate at the level of a factor of ~ 6 due to uncertainty in the initial 
mass function (Kroupa, 2002) and by up to a factor of 2-3 due to uncertainty 
in the other parameters (de Mink and Belczynski, 2015; Klencki et al., 2018). 

The treatment of mass transfer represents a key source of uncertainty in 
isolated binary evolution (Belczynski et al., 2021). Is the mass transfer dynam- 
ically stable? Stability thresholds can be formulated in terms of the mass ratio 
(Claeys et al., 2014) or the relative radial response of the star and its Roche 
lobe to mass transfer (Soberman et al., 1997; Ge et al., 2015). The depth of the 
convective envelope of the donor may play a critical role (Klencki et al., 2021). 


22 Ilya Mandel, Floor S. Broekgaarden 


If the mass transfer is dynamically stable, how conservative is it, i.e., what frac- 
tion of mass lost by the donor does the companion accrete (e.g., Kippenhahn 
and Meyer-Hofmeister, 1977), including the case of accretion onto compact 
objects (e.g., van Son et al, 2020)? And how much angular momentum is 
carried away by any mass that leaves the system, changing the orbital sepa- 
ration of the binary (e.g., Vinciguerra et al., 2020)? How does mass transfer 
in eccentric binaries affect their orbit (e.g.,. Sepinsky et al., 2010; Dosopoulou 
and Kalogera, 2016)? If the mass transfer is dynamically unstable, leading to 
the formation of a common envelope, can the envelope be ejected, and what 
sets the orbital separation after the envelope ejection? Simple energy scalings 
are typically used to treat the common-envelope phase of binary evolution 
(Webbink, 1984; de Kool, 1990), partly because, despite decades of modelling 
common envelopes (Ivanova et al., 2013), the first incomplete models of this 
phase for massive stars are only just appearing (Ricker et al., 2019; Law-Smith 
et al., 2020; Lau et al., 2022). Last but not least, how does the history of mass 
transfer impact the final outcomes of stellar evolution, including supernova 
explosions (e.g., Brown et al., 2001; Schneider et al., 2021)? 

A range of uncertainties in single stellar evolution play a key role. How 
much massive stars expand at various stages in their lives determines the onset 
and outcomes of mass transfer (e.g., Laplace et al., 2020). Stellar winds affect 
the final masses and separations (Vink, 2017). The efficiency of tidal coupling 
determines orbital circularisation (Zahn, 1977; Hut, 1981; Vick and Lai, 2018). 
The amount of mass lost in supernovae and the natal kicks that supernova 
remnants receive impact binary survival during the stellar explosions (Sipior 
and Sigurdsson, 2002; Fryer et al., 2012; Müller, 2020). Even relatively rare 
supernova variants such as electron-capture supernovae and (pulsational) pair- 
instability supernovae are important for coalescing compact object formation 
(e.g., Belezynski et al., 2018b; Stevenson et al., 2019; Farmer et al., 2020). 


4.1.2 Isolated binaries: chemically homogeneous evolution and population III 
stars 


One way to avoid the problem of fitting very extended stars into a binary 
tight enough to merge through gravitational waves is to prevent the stars 
from expanding in the first place. Here we consider two magic wands that 
might prevent expansion: efficient mixing or initially metal-free composition. 

Suppose the stars are sufficiently massive, their metallicity is low enough 
to suppress strong winds, and they are close enough to be tidally locked in an 
orbit that enforces rotation at a significant fraction of the break-up velocity. 
Then the formation of temperature gradients between the poles and equator 
of the star may lead to large-scale circulation and efficient mixing (Eddington, 
1925; Sweet, 1950; Endal and Sofia, 1978). The mixing of helium out of the 
core into the envelope and fresh hydrogen from the envelope into the core 
allows the star to evolve chemically homogeneously, fusing almost all of the 
hydrogen into helium on the main sequence (Heger et al., 2000; Maeder and 
Meynet, 2000; Yoon et al., 2006). After the end of hydrogen fusion, the entire 
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star contracts into a naked helium star. Such systems avoid mass transfer 
after the main sequence, and may form heavy, ultimately coalescing black hole 
binaries in situ (Mandel and de Mink, 2016; Marchant et al., 2016; de Mink and 
Mandel, 2016). The high mass required for chemically homogeneous evolution 
makes this channel relevant only for the formation of merging black holes and, 
perhaps, rare NS-BH binaries (Marchant et al., 2017). 

We can write a similar Drake's equation to Eq. (2) for the chemically homo- 
geneous evolution channel. Following Mandel and de Mink (2016), the fraction 
of binaries that form coalescing binary black holes through this channel is 


fpu-BH = fprimary x fsecondary x fiuit sep X fz x Jerga 
~ 0.0002 x 0.5 x 0.1 x 0.1 x 1 = 1076. (3) 


Only the most massive stars with initial masses above ~ 50 Mọ are likely 
to experience sufficient mixing to evolve chemically homogeneously, lower- 
ing fprimary to 7 2 X 10-4. A flat mass ratio distribution would again yield 
fsecondary © 0.5. The range of initial separations allowing for chemically ho- 
mogeneous evolution is at most a factor of two: stars that are initially too 
far apart either avoid tidal locking or do not rotate sufficiently rapidly, while 
stars that are too close overflow the L2 Lagrange point and promptly merge. If 
initial separations cover some 5 orders of magnitude and a uniform in the log- 
arithm distribution of separations is assumed, up to finit sep © 0.1 of binaries 
with the right mass range could still experience chemically homogeneous evo- 
lution. However, this also requires a sufficiently low metallicity, in part to avoid 
binary widening through excessive mass loss in winds, necessitating a factor 
describing the fraction of star formation at metallicities of interest, estimated 
here as fz ~ 0.1. On the other hand, if a binary evolves chemically homoge- 
neously and forms black holes in situ, it is likely to be sufficiently compact 
to merge within the age of the Universe through gravitational-wave emission, 
SO fmerge ^ 1. The overall yield estimate of fBH-BH 7 10-9 would suggest a 
binary black hole coalescence rate of ~ 10 Gpc^? yr! for this channel. 

The chemically homogeneous evolution channel as described above does 
not involve mass transfer after the main sequence and hence avoids some of 
uncertainties discussed in the previous subsection. On the other hand, the very 
possibility of efficient chemical mixing is uncertain. The observational evidence 
ranges from tentative to contradictory (Almeida et al., 2015; Mandel and de 
Mink, 2016; Abdul-Masih et al., 2021) and theoretical models of mixing differ 
in the predicted thresholds on its efficiency. Moreover, much of the parameter 
space in this channel likely requires over-contact binaries (Marchant et al., 
2016; du Buisson et al., 2020; Riley et al., 2021), which bring extra modelling 
challenges and uncertainties. On the other hand, a star could be spun up 
to rapid rotation and chemically homogeneous evolution by mass accretion 
from the companion (Cantiello et al., 2007). This mode of quasi chemically 
homogeneous evolution could also make a significant contribution to high- 
mass BH-BH mergers originating in low-metallicity environments (Eldridge 
and Stanway, 2016). 
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First-generation metal-free stars in the early Universe, known as popula- 
tion III stars, may represent an alternative way to avoid radial expansion after 
the main sequence. The absence of primordial carbon, nitrogen and oxygen in 
these stars prevents the CNO cycle from operating, so that hydrogen fusion 
can proceed only through the pp-chain reaction. While the core becomes hot 
enough for helium produced in the core to form carbon through the triple- 
alpha reaction, eventually allowing CNO cycle hydrogen fusion to take over, 
this may not happen to the same extent in the hydrogen shell around the he- 
lium core after core hydrogen exhaustion. Therefore, population III stars may 
experience reduced radial expansion relative to more metal-rich stars (Marigo 
et al., 2001) and can avoid mass transfer and many of the complications of 
the classical isolated binary evolution channel. Moreover, this means that the 
hydrogen envelope of such stars remains tightly bound, and is unlikely to be 
ejected during the stellar collapse as would be the case for population I and II 
supergiants (Nadezhin, 1980; Lovegrove and Woosley, 2013; Fernández et al., 
2018) allowing the star to retain most of its mass (Kinugawa et al., 2021), 
akin to chemically homogeneously evolving stars. This mass retention is fur- 
ther aided by the drastically reduced wind-driven mass loss at zero metallicity. 
Finally, the lack of metals has been conjectured to suppress the fragmentation 
of a giant molecular cloud during star formation, increasing the typical mass 
of population III stars and making these prime progenitor candidates for bi- 
nary black hole coalescences (Belczynski et al., 2004; Kinugawa et al., 2014; 
Inayoshi et al., 2016, 2017). 

'The lack of observational evidence and challenges in theoretical models put 
many of these assumptions in doubt. There are conflicting models for the initial 
properties of population III binaries (Hirano et al., 2014; Stacy et al., 2016) 
which could drastically reduce the predicted merger rate (Hartwig et al., 2016; 
Belczynski et al., 2017), while the possibility of significant stellar expansion 
could remove this altogether as an independent channel. 


4.2 Dynamical formation 


Dynamical formation broadly encompasses compact object mergers in which 
gravitational interactions with other stars played a significant role in bringing 
the two merging objects closer. Dynamical formation can be broadly divided 
into several categories based on environment. Historically, the focus has been 
primarily on young star clusters and globular clusters, where dynamical for- 
mation has been explored by Sigurdsson and Hernquist (1993); Kulkarni et al. 
(1993); Portegies Zwart and McMillan (2000); O'Leary et al. (2006); Banerjee 
et al. (2010); Downing et al. (2011); Morscher et al. (2015) through a mix of 
numerical gravitational experiments and semi-analytical estimates. We illus- 
trate the expectations for globular cluster dynamical formation with another 
Drake's equation before discussing other dynamical formation environments: 
galactic nuclei and hierarchical triple systems. 
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4.2.1 Young star clusters and globular clusters 


The compact objects in the cluster may or may not have started out in the 
same binary. In the classical version of dynamical formation in a dense stellar 
environment the black holes, being heavier than other cluster objects, sink 
toward the centre of the cluster through mass segregation (Spitzer, 1969; Bin- 
ney and Tremaine, 2008). Once there, they readily form binaries or substitute 
into existing binaries through three-body interactions: in a binary-single inter- 
action, the lightest object is most likely to be ejected, while the two heavier 
objects remain behind in a binary (Sigurdsson and Phinney, 1993; Ginat and 
Perets, 2021). If the binary is hard — i.e., if its orbital velocity is larger than 
the velocity dispersion for the case of comparable component and interloper 
masses — subsequent interactions will further harden the binary as the energy 
of the binary is shared with the interlopers. Eventually, the binary may harden 
enough for gravitational waves to take over, provided the interaction rate is 
sufficiently high and the binary is not prematurely ejected from the cluster. 
'This suggests that the number of BH-BH mergers per globular cluster is: 


NBH-BH, cluster — Nstars x feu x fretain x fmerge 
~ 109 x 0.001 x 0.2 x 1 = 200. (4) 


Here, we focus on BH-BHs, because neutron stars, being lighter, are pref- 
erentially ejected during interactions, reducing their formation rate (Ye et al., 
2020). We consider the heaviest globular clusters with Netars © 10° stars. We 
optimistically assume that all > 20Mo stars will form a BH, ignoring the 
likelihood that natal kicks will eject some of the BHs from the cluster, whose 
escape velocity is only ~ 50km s-!. This sets fgg œ% 0.001. In order for a 
binary to merge through the emission of gravitational-waves within 10 Gyr, its 
orbital velocity must be ~ 500 km s^! (Peters, 1964). This is almost two orders 
of magnitude larger than the typical velocity dispersion of 10 km s~!. Tens of 
strong binary-single scattering interactions, mostly with other black holes, will 
be necessary to harden one binary by two orders of magnitude (Quinlan, 1996). 
Once the binary becomes sufficiently hard, the lightest of the three interact- 
ing objects will typically be ejected from the cluster, depleting the fraction of 
retained black holes by a factor of a few (Heggie and Hut, 2003). We therefore 
set the fraction of black holes retained for possible mergers to fretain 7 0.2. Fi- 
nally, we consider a cluster sufficiently dense (central BH density above ~ 10° 
pc ?) that most of the BHs that avoid being ejected do succeed in merging 
within 10 Gyr, so that fmerge = 1. Even with these optimistic assumptions, 
and allowing for a constant merger rate over the 10 Gyr lifetime of globular 
clusters, whose space density is ~ 1 Mpc^?, we estimate a maximum merger 
rate of 200 x 10? Gpc ? x 10719 yr-! = 20 Gpc ? yrat. 

'This very crude estimate helps highlight some of the issues that may im- 
pact the event rates. What is the fraction of compact objects retained in 
clusters despite natal kicks? What were the initial properties (mass distribu- 
tion, binary fraction, density, velocity dispersions) of globular clusters, before 
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the = 10 Gyr of evolution typical for the locally observed clusters (Fragione 
and Kocsis, 2018)? What fraction of merger products avoid being ejected by 
gravitational-wave recoil kicks and can be reused for subsequent hierarchical 
mergers (Rodriguez et al., 2018; Gerosa and Fishbach, 2021)? What role does 
ongoing stellar evolution, including BH growth through tidal disruption of 
stars and gas accretion (Vesperini et al., 2010), play in dynamical formation? 
While multiple approaches to modelling cluster gravitational dynamics (e.g., 
N-body and Monte Carlo models) are expected to yield consistent results, the 
range of predictions highlighted in the executive summary tables and figures 
illustrates the impact of different approaches to answering the questions posed 
in this paragraph. 


4.2.2 Nuclear star clusters 


Dynamical formation has also been explored in the context of nuclear clusters 
in the centres of galaxies. Miller and Lauburg (2009) and Antonini and Rasio 
(2016) pointed out that nuclear clusters in the absence of a massive black hole, 
which are characterised by deeper potentials and hence higher escape velocities 
than globular clusters, could be promising environments for BH-BH dynamical 
formation. O'Leary et al. (2009); Tsang (2013); Hoang et al. (2020); Rasskazov 
and Kocsis (2019) explored direct gravitational-wave bremsstrahlung captures 
from 2-body scattering. Other potential contributions include the deposition 
of stellar-mass black holes near a massive black hole by the infall of a globular 
cluster (Arca-Sedda and Gualandris, 2018). We quote a selection of the latest 
predictions for all of these channels in the executive summary. 

Meanwhile, Bellovary et al. (2016); Bartos et al. (2017); Stone et al. (2017); 
McKernan et al. (2018) emphasised the important role that an accretion disk 
in an active galactic nucleus could play in enhancing the rate of binary black 
hole mergers. The rich physics of the deep potential well combined with gas 
interactions can yield new merger channels (Tagawa et al., 2020). Compact ob- 
jects may be captured by the disk or directly form in the disk; in fact, gas may 
not even be strictly necessary for this, as vector resonant relaxation against 
the stellar background could naturally lead massive objects to settle into a dis 
(Szólgyén and Kocsis, 2018). Migration traps have been suggested as breed- 
ing grounds for dynamical interactions that will enhance the BH+BH merger 
rate (Bellovary et al., 2016; Yang et al., 2019; Secunda et al., 2019; Peng and 
Chen, 2021), although there is robust debate in the literature about their ef- 
ficacy (e.g., Pan and Yang, 2021). Gas accretion may aid gravitational waves 
in driving binaries toward merger as well as grow the black hole mass. Mean- 
while, the very high escape velocity in the vicinity of a supermassive black hole 
ensures that recoil kicks are unlikely to eject merger products, so hierarchical 
mergers should be more common than in other dynamical environments (Yang 
et al., 2019; Secunda et al., 2020; Tagawa et al., 2021). The rich physics of ac- 
tive galactic nuclei make this channel both promising and complex, with many 
details such as the feedback from BH-gas interaction and accretion remaining 
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a topic of ongoing work, leading to a broad range of merger rate predictions 
(Tagawa et al., 2020). 


4.2.8 Hierarchical 3-body systems 


Another notable dynamical formation channel involves hierarchical triple sys- 
tems (see Naoz 2016 for a review). In these systems, the transfer of angular 
momentum back and forth between the inner binary and the wider outer orbit 
of the third companion can cyclically increase the eccentricity of the inner 
binary (Lidov, 1962; Kozai, 1962). Once the eccentricity of the inner binary 
becomes sufficiently high, energy can be dissipated by gravitational waves or 
tides. This could drive the inner binary to merger (e.g., Wen, 2003; Liu and 
Lai, 2018). This channel may operate in a number of regimes, from isolated 
stellar triples (Silsbee and Tremaine, 2017) to dynamically formed triples in 
clusters (Antonini et al., 2016) to triples in which the outer companion is a 
massive black hole in a galactic center (Antonini and Perets, 2012; Hoang 
et al., 2018). 

Unlike other dynamical formation scenarios, which are most promising for 
BH-BH mergers, hierarchical triples could be important for NS-NS forma- 
tion, although predictions are sensitive to assumed NS natal kicks (Hamers 
and Thompson, 2019). Extensions to the triple channel include higher mul- 
tiplicity systems, such as quadruples (Fragione et al., 2020; Hamers et al., 
2021; Vynatheya and Hamers, 2021). Some authors have associated specific 
gravitational-wave sources, such as GW190814, with formation in hierarchical 
triples (e.g., Lu et al., 2021). 

The triple channel can also operate in conjunction with other channels. For 
example, stellar evolution and associated mass loss could force a previously sta- 
ble stellar triple into instability or drive the inner binary toward merger (Perets 
and Kratter, 2012; Shappee and Thompson, 2013; Stephan et al., 2016), while 
chemically homogenous evolution and triple dynamics could produce sequen- 
tial BH-BH mergers (Vigna-Gómez et al., 2021b). In fact, it is the combination 
of triple dynamics with stellar evolution, mass transfer, tides and gravitational 
waves (as well as the possibility of other dynamical interactions with external 
perturbers) which makes it particularly challenging to analyse. The uncertain 
initial conditions, such as a possible correlation between the mutual inclina- 
tion, separations and eccentricities of the inner and outer binaries (Moe and 
Di Stefano, 2017), query the efficacy of this channel (Toonen et al., 2020) (but 
see Rose et al. 2019, who find that the observed orbits of massive stars may 
be consistent with being driven by triple dynamics). 

Field binaries perturbed by flybys from other stars represent a special case 
on the border between isolated and dynamically formed binaries. Such flybys 
can enhance the eccentricity of a wide compact-object binary that evolved as 
an isolated system, driving it to merge by gravitational-wave emission (Raveh 
et al., 2022; Michaely and Naoz, 2022). 
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4.3 Exotica 


So far, we focussed exclusively on black holes of astrophysical origin, rather 
than primordial black holes forming from the collapse of early-Universe per- 
turbations. The latter are potentially very interesting sources of gravitational 
waves. Bird et al. (2016); Ali-Haimoud et al. (2017); Chen and Huang (2018) 
have argued that a significant fraction of dark matter could be contained in 
black holes with masses of tens of solar masses, and their mergers could be 
responsible for some of the gravitational-wave signals observed to date. How- 
ever, even allowing for a significant amount of mass in primordial black holes, 
merger rate predictions are sensitive assumptions about the initial distribu- 
tion of such sources in binaries or small clusters (e.g., Ali-Haimoud et al., 2017; 
Korol et al., 2020; De Luca et al., 2020a) as well as possible accretion onto 
primordial black holes (De Luca et al., 2020b). Meanwhile, some of the sce- 
narios proposed for the formation of BH-BH binaries involve physics beyond 
the standard model (Sakstein et al., 2020) and cosmological coupling (Croker 
et al., 2021). 

The discovery of the BH-BH merger GW190521, which left behind a ~ 
150 Ms remnant (Abbott et al., 2020c), naturally renews interest in the possi- 
bility of intermediate mass black hole (IMBH) mergers. IMBH mergers may be 
a natural consequence of binary or dynamical evolution (e.g., Amaro-Seoane 
and Santamaría, 2010; McKernan et al., 2012; Belczynski et al., 2014), and 
the gravitational waves from mergers of few hundred solar-mass IMBHs could 
be detected by current instruments (Veitch et al., 2015; Graff et al., 2015; 
Abbott et al., 2017d). Alternatively, stellar-mass compact objects could spiral 
into IMBHs in clusters. These intermediate-mass-ratio inspirals are also po- 
tentially detectable as gravitational-wave signals (Mandel et al., 2008; Haster 
et al., 2016b,a). However, the prevalence of IMBHs remains highly uncertain 
despite a variety of approaches to their detection (e.g., Pasham et al., 2014; 
Paynter et al., 2021). The challenges associated with observing IMBHs are 
discussed by Miller and Colbert (2004) and Greene et al. (2020). Given the 
many orders of magnitude uncertainties, we resist the temptation to review 
the merger rates of binaries involving IMBHs. Similarly, we do not specifically 
discuss hierarchical mergers of stellar-mass black holes or IMBHs as a channel 
for forming supermassive back holes at high redshifts (e.g., Volonteri, 2010; 
Kroupa et al., 2020). The heaviest merging binaries observed through gravita- 
tional waves could represent the tail end of this process of hierarchical massive 
black-hole formation associated with mergers of ultra-dwarf galaxies (Palmese 
and Conselice, 2021), and many more such mergers may be observable with 
third-generation gravitational-wave detectors (Gair et al., 2009). 


5 Outlook 


We reviewed the current estimates of compact-binary merger rates, includ- 
ing both direct observational constraints and theoretical models motivated by 
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observations. We highlighted some of the key sources of uncertainty in both 
observations and models. Of course, these uncertainties can be alternatively 
viewed as opportunities. 

The number of observations is growing rapidly. The number of direct 
gravitational-wave detections has grown from 0 to 50 between 2015 and 2020. 
The range, the number of observed sources, and the measurement accuracy will 
all increase as the detector sensitivity improves (Abbott et al., 2020b). Eventu- 
ally, third-generation gravitational-wave instruments will be able to precisely 
probe the evolving compact binary merger rate and mass distribution across 
cosmic time (e.g., Kalogera et al, 2019; Ng et al., 2021). The number and 
quality of electromagnetic observations are also rapidly growing, again with a 
very positive outlook, e.g., for kilonovae as new instruments like the Vera C. 
Rubin Observatory come online (e.g., Chase et al., 2021). 

These improvements in observations are being matched by a growing so- 
phistication in the models. On the one hand, increasingly complex models of 
binary evolution and dynamics are sensitive to a broad range of uncertain as- 
sumptions about initial conditions, stellar evolution, mass loss, mass transfer, 
supernova kicks, and the cosmic star formation history of the Universe. On 
the other hand, there is growing evidence that observations of the rates and 
properties of compact-object mergers, coupled with other observational con- 
straints, will make it possible to resolve these uncertainties and come away 
with a deep understanding of the physics of compact binary formation (e.g., 
Mandel and O'Shaughnessy, 2010; Stevenson et al., 2015; Barrett et al., 2018; 
Wong et al., 2021). 

While we focussed on overall rate estimates in this review, specific forma- 
tion channels are likely to carry distinguishable observational signatures. We 
mention two predictions for merging binary black holes by way of example. 
Dynamically formed binary black holes should have isotropically oriented spin 
directions, while BH-BH mergers from isolated binary evolution are expected 
to have spins preferentially aligned with the orbital angular momentum (e.g., 
Farr et al, 2017; Rodriguez et al., 2016b). Is there already evidence for a 
separate population of dynamically formed systems based on spin-orbit mis- 
alignment angles (Abbott et al., 2021c; Roulet et al., 2021; Galaudage et al., 
2021)? Meanwhile, the orbits of most binary black holes are expected to circu- 
larise through gravitational-wave emission by the time they approach merger 
and become detectable with current instruments (Peters, 1964). But a ~ 596 
fraction of binaries that evolve through binary-single dynamical scatterings 
may buck the trend, becoming captured at very high eccentricities and small 
separations through gravitational-wave bremsstrahlung during chaotic inter- 
actions (Samsing et al, 2014; Rodriguez et al., 2018). Could the tentative 
evidence for eccentric binary black holes in gravitational-wave data (Romero- 
Shaw et al., 2021) lead to an inferred measurement of the rate of dynamically 
formed binaries (Zevin et al., 2021)? 

Reviewing such a rapidly developing field presents obvious challenges, and 
we are aware that we have only shown a selection of highlights. We hope that 
this is still a useful snapshot of the status of the field against which future 
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rapid progress can be compared. In fact, we find ourselves in the awkward 
position of wishing for this review to become obsolete soon after it is written, 
to be updated as the field moves forward! 
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